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SUMMARY 

A classical view of speciation is that reproductive isola- 
tion arises as a by-product of genetic divergence. Here, 
individual-based simulations are used to evaluate whether 
the mechanisms implied by this view may result in rapid 
speciation if the only source of genetic divergence are mu- 
tation and random genetic drift. Distinctive features of 
the simulations are the consideration of the complete pro- 
cess of speciation (from initiation until completion), and 
of a large number of loci, which was only one order of mag- 
nitude smaller than that of bacteria. It is demonstrated 
that rapid speciation on the time scale of hundreds of gen- 
erations is plausible without the need for extreme founder 
events, complete geographic isolation, the existence of dis- 
tinct adaptive peaks or selection for local adaptation. The 
plausibility of speciation is enhanced by population sub- 
division. Simultaneous emergence of more than two new 
species from a subdivided population is highly probable. 
Numerical examples relevant to the theory of centrifugal 
speciation and to the conjectures about the fate of "ring 
species" and "sexual continuums" are presented. 

1. INTRODUCTION 

In spite of significant advances in both theoretical and 
experimental studies of evolution, understanding specia- 
tion still remains a major challenge (Mayr 1982; Coyne 
1992). In recent years a view of speciation that can be 
traced to the classical Dobzhansky model (Dobzhansky 



1937; MuUer 1942) has become a point of renewed in- 
terest (Nei et al. 1983; Wagner et al. 1994; Orr 1995; 
Orr & Orr 1996; Gavrilets 1997a,b; Gavrilets & Hastings 
1996; Gavrilets & Gravner 1997). Using the metaphor of 
"adaptive landscapes" (Wright 1932), populations diverge 
genetically along the "ridges" of highly-fit genotypes and 
become reproductively isolated species when they come 
to be on opposite sides of a "hole" in a "holey" adaptive 
landscape (Gavrilets 1997a; Gavrilets & Gravner 1997). 

Current renewal of interest in Dobzhansky-type mod- 
els is coming from different sources. First, there is grow- 
ing experimental support for the genetic architecture im- 
plied (reviewed in Orr 1995; Orr & Orr 1996; Gavrilets 
1997a, b). On the theoretical side, there is a growing real- 
ization of serious problems with Sewall Wright's metaphor 
of "rugged" adaptive landscapes (Whitlock et al. 1995; 
Gavrilets 1997a) and shifting balance theory (Gavrilets 
1996; Coyne et al. 1997). Both have for a long time 
have been a focus of theoretical studies of evolution. On 
the other hand, it recently has been demonstrated that 
the existence of "ridges" of highly-fit genotypes, which 
is postulated in Dobzhansky-type models, is actually a 
general property of multidimensional adaptive landscapes 
(Gavrilets & Gravner 1997; Reidys et al. 1997). 

Earlier theoretical studies of speciation as a conse- 
quence of genetic divergence along "ridges" in the adap- 
tive landscape have shown it to be a distinct possibility. 
These studies involved various approximations such as 
small number of loci (Nei et al. 1983; Wagner et al. 1994; 
Gavrilets & Hastings 1996), linkage equilibrium (Nei et 
al. 1983), infinitely large number of highly mutable un- 
linked loci and preferential mating (Higgs & Derrida 1992; 
Manzo & Peliti 1994), or absence of within-population 
variability during divergence (Orr & Orr 1996). Here we 
extend the previous work and report results of individual- 
based simulations that did not involve any approxima- 
tions typical in studying multilocus systems. Using a 
linear algorithm for generating random numbers with a 
given distribution (Vose 1991) has allowed us to model 
hundreds of linked loci potentially influencing fitness. 
The number of considered loci, which is only one order 
of magnitude smaller than that of bacteria (Blattner et 
al. 1997), is a significant distinctive feature of our sim- 
ulations. An analytic theory for the dynamic patterns 
observed will be published elsewhere (Gavrilets unpub- 
lished) . 

2. MODEL 

We model finite subdivided populations of sexual haploid 
individuals different with respect to L diallelic loci. Evo- 
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lutionary factors included were mutation, recombination, 
migration, genetic drift and selection. To decrease the 
number of parameters we used some symmetry assump- 
tions. Mutation rates were equal for forward and back- 
ward mutations and across loci. Recombination rates be- 
tween adjacent loci on the same chromosome were equal 
and recombination events took place independently of 
each other. The population was subdivided into subpop- 
ulations of equal size. 

Without any loss of generality each individual can be 
represented as a sequence of O's and I's. Let 1°' = 
where If is equal to or 1, be such a sequence 
for an individual a. We will be interested in the levels 
of genetic variability within subpopulations and genetic 
divergence between subpopulations. Both can be charac- 
terized in terms of genetic distance d defined as the num- 
ber of loci at which two individuals are different. More 
formally, the genetic distance d"'^ between individuals a 
and /3 is 

^a^^j-ilT-l^r. (1) 

Genetic distance d is the standard Hamming distance. To 

characterize the levels of genetic variability within sub- 
populations and genetic divergence between subpopula- 
tions we used the average genetic distances (AGD) be- 
tween individuals from the same subpopulation and from 
different subpopulations, respectively. To reflect the idea 
that reproductive isolation arises simultaneously with ge- 
netic divergence we posit that an encounter of two in- 
dividuals can result in mating and viable and fecund off- 
spring only if the individuals are different in no more than 
K loci (cf.. Higgs & Derrida 1992). Otherwise the indi- 
viduals do not mate (premating reproductive isolation) 
or their offspring is inviable or sterile (postmating repro- 
ductive isolation) . More formally, we assign "fitness" w to 
each pair of individuals depending on the genetic distance 
d between them 

( 1 {or d<K, , , 

^('^^ = I {0Td>K. 

In this formulation, any two genotypes different in more 
than K loci can be considered as sitting on opposite sides 
of a hole in a holey adaptive landscape. At the same 
time, a population can evolve to any reproductively iso- 
lated state by a chain of single locus substitutions. The 
neutral case (no reproductive isolation) corresponds to K 
equal to the number of loci. In our simulations individ- 
uals migrated before mating and pairing of individuals 
occurred randomly. There were 384 on a single chromo- 



some, the rate of recombination between adjacent loci was 
0.005, and the mutation rate per locus was lO"'^. The K 
values used are 10, 20 and 30 which are within the range 
of estimates of the minimum number of genes involved in 
reproductive isolation (Singh 1990; Coyne and Orr 1998). 

Our simulations started with all N individuals identi- 
cal. During the first 1000 generations there were no re- 
strictions on migration between subpopulations and the 
whole population evolved as a single randomly mating 
unit. One thousand generations was sufficient for the 
population to reach a state of stochastic equilibrium be- 
tween factors affecting its dynamics. This was appar- 
ent from the behavior of AGD (see figures below) as well 
as from the behavior of the number of segregating loci 
and the number of different genotypes (data not pre- 
sented). Starting with generation 1000, restrictions on 
migration wctc introduced. We used both the island 
model (\\ right 1931) and the one-dimensional stepping- 
stone niodcil (Kimura & Weiss 1964). In the former case, 
with probability m an individual migrated to another sub- 
population (with all subpopulations equally accessible). 
In the later ease, the subpopulations were arranged along 
a line with migration only between neighboring subpop- 
ulations (with rate m/2). The simulations were contin- 
ued until generation 3000. We have chosen this relatively 
short time span for two reasons. First, here we are inter- 
ested in the plausibility of rapid speciation. Second, only 
during relatively short time intervals may one assume the 
constancy of abiotic and biotic environment implied in 
our model. 

3. NUMERICAL RESULTS 

Extensive simulations have revealed the existence of three 

different dynamical regimes. In the first regime, the AGD 
both between and within subpopulations reach (stochas- 
tic) equilibrium values below K (Fig. la). Between- 
subpopulations matings take place and are fecund. No 
further genetic divergence and speciation take place (at 
least until generation 3000). In the second regime, the 
AGD between subpopulations steadily increases after in- 
troduction of restrictions on migration (Fig. lb). Subpop- 
ulations diverge genetically accumulating different muta- 
tions. After some time, three effects take place simulta- 
neously: (1) AGD between subpopulations significantly 
exceed AGD within them, (2) encounters between indi- 
viduals from different subpopulations do not result in vi- 
able and fertile offspring, (3) evolutionary changes in a 
subpopulation do not affect other subpopulation. Thus, 
subpopulations form separate genotypic clusters in geno- 
type space, become reproductively isolated and undertake 
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changes as evolutionary independent units. These dy- 
namics are interpreted as speciation according to any of 
the species concept common in the hterature (e.g., Mallet 
1995; Claridge et al. 1997). Here, speciation has occurred 
in the presence of gene flow and, thus, is parapatric. The 
third regime is a combination of the first two. After some 
transient time during which the AGD both within and 
between subpopulations have appeared to reach a new 
stochastic equilibrium, the AGD between subpopulations 
starts suddenly increasing (Fig.lc). The system imder- 
goes a stochastic transition resulting in speciation and 
further accumulation of genetic differences between sub- 
populations. After speciation has taken place, the genetic 
variability within (sub) populations decreases. [Note that 
instead of the three different dynamical regimes one may 
choose to think of a single regime with a varying (from 
zero to infinity) time until the initiation of divergence.] 

Dynamic behavior depends on the intensity of different 
factors operating in the system. As expected, reduced mi- 
gration and increased mutation both increase the plausi- 
bility of speciation (data not presented). Population size 
and the strength of selection play a double role. With 
strong selection (small K), increased population size de- 
lays speciation (compare Fig. 2a and b). With weak se- 
lection (large K), decreasing population size delays spe- 
ciation (compare Fig. 2c and d). With smaller population 
sizes, decreasing the strength of selection delays specia- 
tion (compare Fig. 2a and c). With larger population size, 
increasing the strength of selection seems to delay speci- 
ation (compare Fig. 2b and d). 

These effects can be explained in the following way. A 
necessary condition for genetic divergence and speciation 
is fixation of alternative alleles in different subpopula- 
tions. After the (sub)population becomes polymorphic 
at K loci, new mutations are selected against because in- 
dividuals carrying them have a reduced probability of fe- 
cund mating. Genetic drift operating in small populations 
overcomes the effect of selection and allows for genetic 
divergence. Genetic divergence of very large randomly 
mating subpopulations connected by migration will not 
occur. On the other hand, if (sub)populations are too 
small, they lack genetic variability necessary to initiate 
divergence. Genetic divergence of very small randomly 
mating subpopulations connected by migration will not 
occur either. The dynamics of genetic divergence are not 
neutral. During the phases of stochastic equilibrium the 
AGD both within and between subpopulations are well 
below neutral expectations (which are given in Watterson 
1975; Li 1976; Slatkin 1987; Strobeck 1987). After spe- 
ciation has taken place, the rate of between-population 
divergence is smaller than the neutral expectation (which 



is twice the expected number of mutations per generation, 

V). 

Increased population subdivision promotes speciation. 
For example, in a population with N = 800 subdivided 
into two subpopulations no speciation has been observed 
(Fig. 3a). However, subdividing into four "islands" im- 
mediately initiates divergence and emergence of 4 species 
(Fig. 3b). The migration rate in these simulations was 
small (m = 0.01). Speciation can also occur for higher mi- 
gration rates with appropriate spatial subdivision. Fig. 3c 
shows results for a single run of a stepping-stone model 
with 8 subpopulations and m = 0.06 (that is there are 6 
migrants per subpopulation per generation). Here, after 
some transient time the population splits into 2 species. 
Splits have occurred in 70 runs (out of 92 performed), in 
18 of which the population was split into three species. 
Figure 4a shows that the two peripheral populations have 
reduced probability of separating from the whole system. 
Figure 4b suggests that this is so because peripheral pop- 
ulations on average have reduced genetic variability and 
are more similar to their neighbors relative to other pop- 
ulations. Fig. 4c shows that, as expected, genetic distance 
increases with the geographic distance. The AGD be- 
tween the first subpopulation and subpopulations four 
through eight is above K meaning these pairs of subpopu- 
lations are reproductively isolated. However, they evolve 
together and have a common evolutionary trajectory. 

4. DISCUSSION 

Here, individual-based simulations were used to evaluate 

the plausibility of rapid speciation as a result of genetic 
divergence along the ridges of a holey adaptive landscape. 
Distinctive features of the simulations are the consider- 
ation of the complete process of speciation (from initi- 
ation until completion), and of a large number of loci, 
which was only one order of magnitude smaller than that 
of bacteria. Our numerical results have implications for 
several important evolutionary questions that have been 
extensively discussed in the literature but have not been 
approached using explicit mathematical models. 

1. Rapid speciation is plausible without the need for ex- 
treme founder events, complete geographic isolation, the 
existence of distinct adaptive peaks or selection for lo- 
cal adaptation. In our simulations, speciation was ob- 
served when the migration rate was on the order of sev- 
eral individuals per subpopulation per generation. The 
time scale in the simulations was short, meaning that re- 
strictions on migration should not be long-lasting. Sev- 
eral hundreds generations might be sufficient for a signif- 
icant divergence and evolution of reproductive isolation 
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(Fig. 1,2, 3). Restoring migration after that to higher lev- 
els will not return the system back to its initial state of 
free gene exchange between subpopulations. The rates of 
speciation observed here are much higher than in mod- 
els of speciation via transitions between adaptive peaks 
across valleys of maladaptation (e.g. Barton and Rouhani 
1993) and in earlier models of holey adaptive landscapes 
that considered only few loci (Nci et al. 1983; Wagner 
et al. 1994; Gavrilets & Hasting 1996). In our model, 
random genetic drift and mutations were the only source 
of genetic divergence and speciation. Here we tacitly as- 
sumed that the loci controlling reproductive isolation are 
different from those responsible for local adaptation. If 
this were not the case or if there were close linkage be- 
tween the two types of loci, selection for local adaptation 
might accelerate genetic divergence and speciation (e.g. 
Schhiter 1996). 

2. The plausibility of speciation is enhanced by popula- 
tion subdivision. In our model, increasing the number of 
subpopulations makes speciation more plausible (Fig. 3). 
The number of subpopulations can be interpreted as a 
measure of the geographic range size of the population. 
Thus, our results show that populations with larger ge- 
ographic range sizes have a greater probability of speci- 
ation (cf, Rosenzweig 1995). Our conclusion about the 
effect of population subdivision on the probability of spe- 
ciation in Dobzhansky-type models differs from that of 
Orr and Orr (1996). They argued that the degree of pop- 
ulation subdivision has no effect on the rate of speciation 
if speciation is caused by mutation and random drift. The 
problematic assumption of their analysis is that fixation 
of new mutations is a neutral process. However, the ex- 
istence of holes in the adaptive landscape makes the pro- 
cess of substitution non-neutral and new mutations are 
selected against when rare (see above). Such mutations 
are fixed more easily in smaller subpopulations (Fig. 2). 

3. Simultaneous emergence of several new species 
from a subdivided population is highly probable. Hoelzer 
and Melnick (1994) have emphasized that this possibil- 
ity should be incorporated more explicitly in the con- 
temporary methods for reconstructing phylogcnies. Re- 
cently several cases of rapid speciation have been de- 
scribed (Schluter & McPhail 1992; McCune 1996; John- 
son et al. 1996). The most spectacular one is probably 
the origin of hundreds of species of Lake Victoria cichlids 
in 12,000 years (Johnson et al. 1996). Our simulations 
have provided examples of two, three and four new species 
emerging from a subdivided population during a short 
time interval. Numerical examples with hundreds of new 
species can be constructed by increasing the number of 
subpopulations. Overall our simulation model provides 



a justification for a verbal "micro-allopatric" model of 
speciation suggested for cichlid fishes (e.g., Reinthal and 
Meyer 1997). 

4- Peripatric speciation vs. centrifugal speciation. Un- 
der appropriate conditions, a subdivided population will 
split into two or more distinct species. An interesting 
question concerns the likelihood of the split at different 
points along the geographic range of the species (e.g. Gas- 
ton 1998). Mayr's (1952) theory of peripatric speciation 
argues that peripheral populations are more likely to split 
off. Surprisingly, our simulations of a one-dimensional 
stepping-stone model show that the breaks resulting in 
the splitting of peripheral populations have the lowest 
probability (Fig. 4a). A reason for this is that peripheral 
populations have lower genetic variability and are more 
similar to their neighbors than the "inner" subpopula- 
tions (Fig. 4b). That genetic changes resulting in spe- 
ciation are more probable in the middle of the species 
range was argued by Brown (1957) in his theory of cen- 
trifugal speciation. It would be interesting to see if the 
splitting pattern that we have observed stays the same in 
more complex models of spatial arrangement of subpop- 
ulations. 

5. "Ring species" and the instability of the "sexual con- 
tinuum." About one forth of the runs of the stepping- 
stone model with 8 subpopulations did not result in spe- 
ciation within the time interval studied. In these runs, 
the neighboring subpopulations were not reproductively 
isolated. At the same time, the AGD between peripheral 
populations significantly exceeded the threshold value K 
(Fig. 4c) meaning that if these populations somehow came 
into a direct contact they would be reproductively iso- 
lated. The situation observed here is analogous to that in 
"ring species" which are traditionally considered as exam- 
ples of incipient, but incomplete, speciation (e.g., Mayr 
1942, 1952; Wake 1997). The population as a whole is 
characterized by the existence of reproductively isolated 
groups that exchange genes through a series of interme- 
diate genotypes. Maynard Smith and Szathmary (1995, 
pp. 163-167), who called such a population state "sexual 
continuum" , conjectured that "sexual continuum" is un- 
stable and would break up into distinct "species." Noest 
(1997) has proposed a model in which a uniform "sexual 
continuum" was unstable to a range of perturbation. His 
analysis did not address whether the system will evolve to 
a non-uniform "sexual continuum" or to a discontinuous 
state. Although the "sexual continuum" observed in the 
runs described by Fig. 4c appears to be stable, continuing 
simulations for a longer period of time would almost defi- 
nitely result in speciation in all runs. This provides a jus- 
tification for Maynard Smith and Szathmary's conjecture. 
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Developing an analytic theory for predicting the persis- 
tence time of a "sexual continuum" and "ring species" is 
an important direction for future work. 
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Legends 

Figure 1. Three types of dynamics. The population is 
subdivided into two subpopulations. Of the three curves 
apparent after generation 1000, the upper curve gives the 
AGD between the subpopulations, whereas the two re- 
maining curves give the AGD within each subpopulation. 
The total population size A'', and K value, are given on 
the graphs. The arrows mark generation 1000 when re- 
strictions on migration (m = 0.01) we introduced. 

Figure 2. Effects of population size and the strength of 
selection on evolutionary divergence and speciation. The 
population is subdivided into two subpopulations. AGD 
between the subpopulations are shown (30 runs for each 
parameter configuration), (a) With K = IQ and N = 200, 
subpopulations steadily diverged immediately after gen- 
eration 1000 in all runs. (6) However, with A'^ = 400 the 
time until the initiation of divergence varied and diver- 
gence was not observed in 13 runs. Comparing the slopes 
of the curves representing AGD between subpopulations 
in (a) and (6) shows that increasing population size de- 
creases the rate of divergence as well, (c) With if = 30 
and N = 200 speciation was observed in only 2 runs, (d) 
With N = 400, speciation took place in 18 runs. Other 
parameters are the same as in Fig.l. 

Figure 3. Effects of spatial subdivision {N = 800, K = 
20; number of loci and the rates of recombination and 
mutation are the same as in Fig.l). Figures (a) and (6) 
show AGD within and between subpopulations for a sin- 
gle run. Results for the 29 other runs are very similar. 
(a) Population is subdivided into 2 equal subpopulations 
with migration rate m — 0.01. No speciation is observed 
(ef. Fig. la). (6) Population is subdivided into 4 equal 
subpopulations; island model with m — 0.01. All AGD 
between subpopulations steadily increase after generation 
1000 indicating emergence of 4 species (cf. Fig. lb). Fig- 
ure (c) corresponds to a stepping-stone model with 8 local 
subpopulations with migration rate m/2 = 0.03 between 
neighboring subpopulations. AGD between neighboring 
subpopulations for a single run resulting in the emergence 
of two species are shown. All distances remain below K 
(no reproductive isolation between neighboring subpop- 
ulations) except between subpopulations 4 and 5 which 
starts increasing around generation 2300 (cf. Fig.lc) in- 
dicating the split into two species with 4 subpopulations 
each. 

Figure 4. Splitting patterns in a stepping-stone system 
with 8 local subpopulations. All parameters are the same 
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as in Fig. 3c. (a) The observed distribution of the position 
of the first break using data from 70 runs (out of 92 total 
runs) that resulted in speciation. The abscissa gives the 
pair of subpopulations between which the split occurred. 
(6) Shown are the AGD within the subpopulations (the 
bars above the integers from 1 to 8) and between the 
subpopulations (the bars between the integers) calculated 
over generations 1500 through 3000 in 22 runs that did 
not result in speciation. (c) The AGD between the first 
subpopulation and other 7 subpopulations using the same 
data as in (6). 
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